## Replication code for: Donnaloja, V. and M. Vink (2023). 
# Like parent, like child: how attitudes towards immigrants spill over to the political inclusion of their children. 
# Journal of Ethnic and Migration Studies.

# This script contains the code to reproduce Figures S1, S2 and S9 from the Supplementary Materials
##interaction terms

#clear all from environment
rm(list = ls())

library("rio")
library("cregg")
library ("FindIt")
library ("FindIt")
library(ggplot2)
library(broom)
library(haven) 
library(expss)

#set working directory



# set working directory or create Project
#import data
mydata <- read_dta("new_data.dta")

## factorise attributes and levels for the ones that have ordered levels
##the outcome stays numeric because otherwise causalanova does not work
mydata$sex <- rio::factorize(mydata$sex)
mydata$religion <- rio::factorize(mydata$religion)
mydata$residence <- rio::factorize(mydata$residence)
mydata$language <- rio::factorize(mydata$language)
mydata$friends <- rio::factorize(mydata$friends)
mydata$occupation <- rio::factorize(mydata$occupation)
mydata$education <- rio::factorize(mydata$education)
mydata$team_support <- rio::factorize(mydata$team_support)
mydata$children <- rio::factorize(mydata$children)
mydata$country <- rio::factorize(mydata$country)
mydata$legal_status <- rio::factorize(mydata$legal_status)

#factorise resp characteristics

mydata$region_grouped_it <- rio::factorize(mydata$region_grouped_it)
mydata$age_new <- rio::factorize(mydata$age_new)
mydata$educ_new <- rio::factorize(mydata$educ_new)
mydata$voting_18 <- rio::factorize(mydata$voting_18)
mydata$gender <- rio::factorize(mydata$gender)
mydata$income_hml_it <- rio::factorize(mydata$income_hml_it)
mydata$empl_stat_eu <- rio::factorize(mydata$empl_stat_eu)
mydata$legal_status <- rio::factorize(mydata$legal_status)

## Specify the order of each factor
mydata$language <- factor(mydata$language,ordered=TRUE,levels=c("Limited","Sufficient", "Excellent"))
mydata$friends <- factor(mydata$friends,ordered=TRUE,levels=c("All Italian","Both", "All foreign"))
mydata$team_support <- factor(mydata$team_support,ordered=TRUE,levels=c("Italy only","Origin only", "Both"))
mydata$children <- factor(mydata$children,ordered=TRUE,levels=c("One","Two", "Four"))
mydata$education <- factor(mydata$education,ordered=TRUE,levels=c("Elementary","Highschool", "university"))
mydata$residence <- factor(mydata$residence, levels=c("Three years", "Five years", "Ten years", "Twenty years"), ordered=TRUE)
mydata$occupation <- factor(mydata$occupation, levels=c("Both work", "Only mother works", "Only father works", "Neither works"), ordered=TRUE)
mydata$religion <- factor(mydata$religion,ordered=FALSE,levels=c("Muslim","Catholic", "No religion"))                                                    
mydata$sex <- factor(mydata$sex,ordered=FALSE,levels=c("Male","Female")) 
mydata$country <- factor(mydata$country, levels=c("Argentina", "China", "Romania", "Pakistan", "Senegal"), ordered=TRUE)
mydata$legal_status <- factor(mydata$legal_status,ordered=FALSE,levels=c("Irregular","Regular"))

table(mydata$voting_18)

#label variables
# Label the variables
mydata = apply_labels(mydata,
                    legal_status = "Parents' legal status",
                    country = "Parents' country of origin",
                    sex = "Child gender",
                    religion = "Religion",
                    occupation = "Parents' employment status",
                    residence = "Parents' length of residence",
                    children = "Total number of children",
                    team_support = "Parents' team support",
                    language = "Parents'Italian proficiency",
                    education = "Parents' education",
                    friends = "Family friends origin")

## change ref category, doesn't work
#mydata <- within(mydata, religion <- relevel(religion, ref = 1))

## %in% to select columns
## constraints
mydata$constraint1constrained <- ifelse(
  (mydata$country %in% c("Romania")) & (mydata$legal_status %in% c("Irregular"))
  , TRUE, FALSE)
mydata$constraint1unconstrained <- ifelse(
  mydata$country %in% c("Argentina", "China", "Pakistan", "Senegal")
  , TRUE, FALSE)

#replication of Egami interaction 

# Fig S1: interaction number of children and country of origin
fit1 <- CausalANOVA(formula=outcome ~ sex  + children + religion  + residence + occupation + education + language + friends + team_support  + country + legal_status,
                    int2.formula = ~ children:country,
                    data=mydata, cluster=mydata$caseid, nway=2)

summary(fit1)
results1 <-summary(fit1)[[4]]

limit <- max(abs(results1$AMIE)) * c(-1, 1)

plotS1 <- ggplot(results1, aes(y= Level2, x=Level1)) +
  geom_raster(aes(fill=AMIE)) +
  scale_fill_gradientn(colours=c("violetred", "white", "orange"),
                       breaks=c(-0.03, -0.02, -0.01, 0, 0.01, 0.02, 0.03), limit=limit) +
  labs(y = "Parents' origin country\n", x = "\nNumber of children in household")

plot(plotS1)
ggsave("FigS1_AMIE_origin_religion.png", width = 14, height = 14, units = "cm", dpi=600)

setEPS()
postscript("FigS1_AMIE_origin_religion.eps")
plotS1
dev.off()

#------------------------------------------------------------------------

#Fig S2: interaction number of children and parents' religion
fit2 <- CausalANOVA(formula=outcome ~ sex  + children + religion  + residence + occupation + education + language + friends + team_support  + country + legal_status,
                    int2.formula = ~ religion:children,
                    data=mydata, cluster=mydata$caseid, nway=2)
summary(fit2)
results <-summary(fit2)[[4]]

limit <- max(abs(results$AMIE)) * c(-1, 1)

plotS2 <- ggplot(results, aes(y= Level2, x=Level1)) +
  geom_raster(aes(fill=AMIE)) +
  scale_fill_gradientn(colours=c("violetred", "white", "orange"),
                       breaks=c(-0.03, -0.02, -0.01, 0, 0.01, 0.02, 0.03), limit=limit) +
        labs(y = "Parents' religion\n", x = "\nNumber of children")

plot(plotS2)
ggsave("FigS2_AMIE_children_religion.png", width = 14, height = 14, units = "cm", dpi=600)

setEPS()
postscript("FigS2_AMIE_children_religion.eps")
plotS2
dev.off()

#------------------------------------------------------------------------
# Fig S10: interaction religion and occupation - among right wing voters
# select on right-wing voters
right_data <- mydata[ which(mydata$voting_18=='Right-wing party'), ]

fit3 <- CausalANOVA(formula=outcome ~ sex  + children + religion  + residence + occupation + education + language + friends + team_support  + country + legal_status,
                     int2.formula = ~ religion:occupation,
                     data=right_data, cluster=right_data$caseid, nway=2)
summary(fit3)

results3 <-summary(fit3)[[4]]

limit <- max(abs(results3$AMIE)) * c(-1, 1)

plotS9 <- ggplot(results3, aes(y= Level2, x=Level1)) +
  geom_raster(aes(fill=AMIE)) +
  scale_fill_gradientn(colours=c("violetred", "white", "orange"),
                       breaks=c(-0.03, -0.02, -0.01, 0, 0.01, 0.02, 0.03), limit=limit) +
  labs(x = "\nParents' religion", y = "Parents' occupation\n")

plot(plotS10)
ggsave("FigS10_AMIE_religion_occup.png", width = 14, height = 14, units = "cm", dpi=600)

setEPS()
postscript("FigS10_AMIE_religion_occup.eps")
plotS9
dev.off()

###END####